3 Equivalent Water Thickness/Canopy Water Content from Imaging Spectroscopy Data

Summary

In this notebook we will explore how Equivalent Water Thickness (EWT) or Canopy Water Content (CWC) can be calculated from the Earth Surface Mineral Dust Source Investigation (EMIT) L2A Reflectance Product, then we will apply this knowledge to calculate CWC over the Jack and Laura Dangermond Preserve located near Santa Barbara, CA.

Background

Equivalent Water Thickness (EWT) is the predicted thickness or absorption path length in centimeters (cm) of water that would be required to yield an observed spectra. In the context of vegetation, this is equivalent to canopy water content (CWC) in g/cm^2 because a cm^3 of water has a mass of 1g.

CWC can be derived from surface reflectance spectra because they provide information about the composition of the target, including water content. Reflectance is the fraction of incoming solar radiation reflected by Earth’s surface. Different materials reflect varying proportions of radiation based upon their chemical composition and physical properties, giving materials their own unique spectral signature or fingerprint. In particular, liquid water causes characteristic absorption features to appear in the near-infrared wavelengths of the solar spectrum, which enables an estimation of its content.

CWC correlates with vegetation type and health, as well as wildfire risk. The methods used here to calculate CWC are based on the ISOFIT python package. The Beer-Lambert physical model used to calculate CWC is described in Green et al. (2006) and Bohn et al. (2020). It uses wavelength-dependent absorption coefficients of liquid water to determine the absorption path length as a function of absorption feature depth. Of note, this model does not account for multiple scattering effects within the canopy and may result in oversestimation of CWC (Bohn et al., 2020).

The Jack and Laura Dangermond Preserve and its surrounding lands are one one of the last “wild coastal” regions in Southern California. The preserve is over 24,000 acres and consists of several ecosystem types and is home to over 600 plant species and over 200 wildlife species.

More about the EMIT mission and EMIT products.

References

Requirements - NASA Earthdata Account
- No Python setup requirements if connected to the workshop cloud instance!
- Local Only - Set up Python Environment. See setup_instructions.md in the /setup/ folder
- Local Only - Downloaded necessary files. This is done at the end of the 01_Finding_Concurrent_Data notebook.

Learning Objectives
- Calculate CWC of a single pixel
- Calculate CWC of an ROI

Tutorial Outline

3.1 Setup
3.2 Opening EMIT Data
3.3 Extracting Reflectance of a Pixel
3.4 Calculating CWC
3.4.1 Single Point
3.4.2 DataFrame of Points
3.5 Applying Inversion in Parallel Across an ROI

# Install ray on 2i2c
!pip install "ray[default]"
Requirement already satisfied: ray[default] in /srv/conda/envs/notebook/lib/python3.10/site-packages (2.8.1)
Requirement already satisfied: click>=7.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (8.1.7)
Requirement already satisfied: filelock in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (3.12.4)
Requirement already satisfied: jsonschema in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (4.19.1)
Requirement already satisfied: msgpack<2.0.0,>=1.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.0.5)
Requirement already satisfied: packaging in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (23.1)
Requirement already satisfied: protobuf!=3.19.5,>=3.15.3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (4.25.1)
Requirement already satisfied: pyyaml in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (6.0.1)
Requirement already satisfied: aiosignal in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.3.1)
Requirement already satisfied: frozenlist in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.4.0)
Requirement already satisfied: requests in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (2.31.0)
Requirement already satisfied: numpy>=1.19.3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.24.4)
Requirement already satisfied: aiohttp>=3.7 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (3.8.5)
Requirement already satisfied: aiohttp-cors in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.7.0)
Requirement already satisfied: colorful in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.5.5)
Requirement already satisfied: py-spy>=0.2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.3.14)
Requirement already satisfied: gpustat>=1.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.1.1)
Requirement already satisfied: opencensus in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.11.3)
Requirement already satisfied: pydantic<2 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.10.13)
Requirement already satisfied: prometheus-client>=0.7.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.17.1)
Requirement already satisfied: smart-open in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (6.4.0)
Requirement already satisfied: virtualenv<20.21.1,>=20.0.24 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (20.21.0)
Requirement already satisfied: grpcio>=1.42.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.59.3)
Requirement already satisfied: attrs>=17.3.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (23.1.0)
Requirement already satisfied: charset-normalizer<4.0,>=2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (3.2.0)
Requirement already satisfied: multidict<7.0,>=4.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (6.0.4)
Requirement already satisfied: async-timeout<5.0,>=4.0.0a3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (4.0.3)
Requirement already satisfied: yarl<2.0,>=1.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (1.9.2)
Requirement already satisfied: nvidia-ml-py>=11.450.129 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from gpustat>=1.0.0->ray[default]) (12.535.133)
Requirement already satisfied: psutil>=5.6.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from gpustat>=1.0.0->ray[default]) (5.9.5)
Requirement already satisfied: blessed>=1.17.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from gpustat>=1.0.0->ray[default]) (1.20.0)
Requirement already satisfied: typing-extensions>=4.2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pydantic<2->ray[default]) (4.8.0)
Requirement already satisfied: distlib<1,>=0.3.6 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from virtualenv<20.21.1,>=20.0.24->ray[default]) (0.3.7)
Requirement already satisfied: platformdirs<4,>=2.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from virtualenv<20.21.1,>=20.0.24->ray[default]) (3.10.0)
Requirement already satisfied: jsonschema-specifications>=2023.03.6 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from jsonschema->ray[default]) (2023.7.1)
Requirement already satisfied: referencing>=0.28.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from jsonschema->ray[default]) (0.30.2)
Requirement already satisfied: rpds-py>=0.7.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from jsonschema->ray[default]) (0.10.3)
Requirement already satisfied: opencensus-context>=0.1.3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from opencensus->ray[default]) (0.1.3)
Requirement already satisfied: google-api-core<3.0.0,>=1.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from opencensus->ray[default]) (2.14.0)
Requirement already satisfied: idna<4,>=2.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->ray[default]) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->ray[default]) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->ray[default]) (2023.7.22)
Requirement already satisfied: wcwidth>=0.1.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from blessed>=1.17.1->gpustat>=1.0.0->ray[default]) (0.2.6)
Requirement already satisfied: six>=1.9.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from blessed>=1.17.1->gpustat>=1.0.0->ray[default]) (1.16.0)
Requirement already satisfied: googleapis-common-protos<2.0.dev0,>=1.56.2 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (1.61.0)
Requirement already satisfied: google-auth<3.0.dev0,>=2.14.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (2.25.1)
Requirement already satisfied: cachetools<6.0,>=2.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (5.3.2)
Requirement already satisfied: pyasn1-modules>=0.2.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (0.3.0)
Requirement already satisfied: rsa<5,>=3.1.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (4.9)
Requirement already satisfied: pyasn1<0.6.0,>=0.4.6 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pyasn1-modules>=0.2.1->google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (0.5.1)

3.1 Setup

Import the required Python libraries.

# Import Packages
import os
import glob
import earthaccess
import math
import numpy as np
import xarray as xr
from osgeo import gdal
import rasterio as rio
import rioxarray as rxr
from matplotlib import pyplot as plt
import hvplot.xarray
import hvplot.pandas
import holoviews as hv
import pandas as pd
import geopandas as gp
import sys
from modules.emit_tools import emit_xarray, ortho_xr
from modules.ewt_calc import calc_ewt
from scipy.optimize import least_squares

3.2 Open an EMIT File

EMIT L2A Reflectance Data are distributed in a non-orthocorrected spatially raw NetCDF4 (.nc) format consisting of the data and its associated metadata. To work with this data, we will use the emit_xarray function from the emit_tools.py module included in the repository.

Set a filepath for an EMIT L2A reflectance file and open using the emit_xarray function.

fp = '../../shared/2023-VITALS-Workshop-AGU/data/EMIT_L2A_RFL_001_20230401T203751_2309114_002.nc'
ds = emit_xarray(fp,ortho=True)

Similarly to what we did in the previous notebook, we can plot a single band to get an idea of where our scene is. This is the same scene we processed earlier.

emit_layer = ds.sel(wavelengths=850,method='nearest')
emit_layer.hvplot.image(cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(
    title=f"{emit_layer.wavelengths:.3f} {emit_layer.wavelengths.units}", xlabel='Longitude',ylabel='Latitude')

3.3 Extracting Reflectance of a Single Pixel

We can mask out the -.01 values used to represent the region of the spectra with strong atmospheric water vapor absorption features.

ds['reflectance'].data[:,:,ds['good_wavelengths'].data==0] = np.nan

Retrieve the spectra from a sample point by providing a latitude and longitude along with a method using the sel function. This will select the pixel closest to the provided coordinates.

point = ds.sel(latitude=34.5399,longitude=-120.3529, method='nearest')
point
<xarray.Dataset>
Dimensions:           (wavelengths: 285)
Coordinates:
  * wavelengths       (wavelengths) float32 381.0 388.4 ... 2.486e+03 2.493e+03
    fwhm              (wavelengths) float32 ...
    good_wavelengths  (wavelengths) float32 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0
    latitude          float64 34.54
    longitude         float64 -120.4
    elev              float32 130.3
    spatial_ref       int64 0
Data variables:
    reflectance       (wavelengths) float32 0.01718 0.01754 ... 0.02635 0.02455
Attributes: (12/40)
    ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
    summary:                           The Earth Surface Mineral Dust Source ...
    keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
    Conventions:                       CF-1.63
    sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
    instrument:                        EMIT
    ...                                ...
    spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
    geotransform:                      [-1.20992414e+02  5.42232520e-04 -0.00...
    day_night_flag:                    Day
    title:                             EMIT L2A Estimated Surface Reflectance...
    granule_id:                        EMIT_L2A_RFL_001_20230401T203751_23091...
    Orthorectified:                    True

We can plot this information to see the spectra.

point.hvplot.line(x='wavelengths',y='reflectance',color='black').opts(title=f"Latitude: {point.latitude.values:.3f} Longitude: {point.longitude.values:.3f}")

3.4 Calculating CWC

As mentioned in the background we use the surface reflectance to estimate CWC. The unique spectral signatures allow identification and quantification based upon the wavelength-dependent absorption coefficients of liquid water. The EMIT mission has applied similar approaches to identify dust source minerals as well as methane point source emissions. The path length of liquid water absorption can be estimated by utilizing a least squares inversion to minimize the residuals between the EMIT reflectance and the Beer-Lambert model (Green et al.,2006), which relates the wavelength-dependent absorption to the path length the photons are traveling through the material. During the inversion, the path lengths are iteratively adjusted to match the modeled spectra to the EMIT reflectance within the water absorption feature region from 850 to 1100 nm.

First, define a Beer-Lambert Model function that returns the vector of residuals between measured and modeled surface reflectance.

# https://github.com/isofit/isofit/blob/main/isofit/inversion/inverse_simple.py#L514C1-L532C17
def beer_lambert_model(x, y, wl, alpha_lw):
    """Function, which computes the vector of residuals between measured and modeled surface reflectance optimizing
    for path length of surface liquid water based on the Beer-Lambert attenuation law.

    Args:
        x:        state vector (liquid water path length, intercept, slope)
        y:        measurement (surface reflectance spectrum)
        wl:       instrument wavelengths
        alpha_lw: wavelength dependent absorption coefficients of liquid water

    Returns:
        resid: residual between modeled and measured surface reflectance
    """

    attenuation = np.exp(-x[0] * 1e7 * alpha_lw)
    rho = (x[1] + x[2] * wl) * attenuation
    resid = rho - y

    return resid

We need some lab measurements of the complex refractive index of liquid water to obtain the wavelength-dependent absorption coefficients. They are calculated by taking four times the product of Pi and the imaginary part of the refractive index, divided by wavelength. The refractive index of liquid water per wavelength is provided by the k_liquid_water_ice.csv in the data folder. We can also preview this data to get a better understanding of the information we are using.

wp_fp = '../data/k_liquid_water_ice.csv'
k_wi = pd.read_csv(wp_fp)
k_wi.head()
wvl_1 T = 22°C wvl_2 T = -8°C wvl_3 T = -25°C wvl_4 T = -7°C wvl_5 T = 25°C (H) wvl_6 T = 20°C wvl_7 T = 25°C (S) Index
0 666.7 2.470000e-08 NaN NaN NaN NaN 660.0 1.660000e-08 650.0 1.640000e-08 650.0 1.870000e-08 650.12971 1.674130e-08 0
1 667.6 2.480000e-08 NaN NaN NaN NaN 670.0 1.890000e-08 675.0 2.230000e-08 651.0 1.890000e-08 654.63616 1.777420e-08 1
2 668.4 2.480000e-08 NaN NaN NaN NaN 680.0 2.090000e-08 700.0 3.350000e-08 652.0 1.910000e-08 660.69347 1.939950e-08 2
3 669.3 2.520000e-08 NaN NaN NaN NaN 690.0 2.400000e-08 725.0 9.150000e-08 653.0 1.940000e-08 665.27314 2.031380e-08 3
4 670.2 2.530000e-08 NaN NaN NaN NaN 700.0 2.900000e-08 750.0 1.560000e-07 654.0 1.970000e-08 669.88461 2.097930e-08 4
fig, axs = plt.subplots(2,4, figsize=(15, 6),  sharex=True, sharey=True, constrained_layout=True)
axs = axs.ravel()
col_n = 0
for i in range(0, 7):
    x = k_wi.iloc[:, col_n+i]
    y = k_wi.iloc[:, col_n+i+1]
    axs[i].scatter(x, y)
    axs[i].set_title(y.name)
    col_n+=1
fig.supylabel('imaginary parts of refractive index')
fig.supxlabel('wavelength')
plt.show()

Now define a function to get the desired data from the csv file.

# https://github.com/isofit/isofit/blob/dev/isofit/core/common.py#L461C1-L488C26
def get_refractive_index(k_wi, a, b, col_wvl, col_k):
    """Convert refractive index table entries to numpy array.

    Args:
        k_wi:    variable
        a:       start line
        b:       end line
        col_wvl: wavelength column in pandas table
        col_k:   k column in pandas table

    Returns:
        wvl_arr: array of wavelengths
        k_arr:   array of imaginary parts of refractive index
    """

    wvl_ = []
    k_ = []

    for ii in range(a, b):
        wvl = k_wi.at[ii, col_wvl]
        k = k_wi.at[ii, col_k]
        wvl_.append(wvl)
        k_.append(k)

    wvl_arr = np.asarray(wvl_)
    k_arr = np.asarray(k_)

    return wvl_arr, k_arr

Lastly, to calculate CWC we define a function that uses least squares optimization to minimize the residuals of our Beer-Lambert Model and find a likely path length of liquid water.

# https://github.com/isofit/isofit/blob/main/isofit/inversion/inverse_simple.py#L443C1-L511C24
def invert_liquid_water(
    rfl_meas: np.array,
    wl: np.array,
    l_shoulder: float = 850,
    r_shoulder: float = 1100,
    lw_init: tuple = (0.02, 0.3, 0.0002),
    lw_bounds: tuple = ([0, 0.5], [0, 1.0], [-0.0004, 0.0004]),
    ewt_detection_limit: float = 0.5,
    return_abs_co: bool = False,
):
    """Given a reflectance estimate, fit a state vector including liquid water path length
    based on a simple Beer-Lambert surface model.

    Args:
        rfl_meas:            surface reflectance spectrum
        wl:                  instrument wavelengths, must be same size as rfl_meas
        l_shoulder:          wavelength of left absorption feature shoulder
        r_shoulder:          wavelength of right absorption feature shoulder
        lw_init:             initial guess for liquid water path length, intercept, and slope
        lw_bounds:           lower and upper bounds for liquid water path length, intercept, and slope
        ewt_detection_limit: upper detection limit for ewt
        return_abs_co:       if True, returns absorption coefficients of liquid water

    Returns:
        solution: estimated liquid water path length, intercept, and slope based on a given surface reflectance
    """
    
    # Ensure least squares is done with float64 datatype
    wl = np.float64(wl)
    
    # params needed for liquid water fitting
    lw_feature_left = np.argmin(abs(l_shoulder - wl))
    lw_feature_right = np.argmin(abs(r_shoulder - wl))
    wl_sel = wl[lw_feature_left : lw_feature_right + 1]

    # adjust upper detection limit for ewt if specified
    if ewt_detection_limit != 0.5:
        lw_bounds[0][1] = ewt_detection_limit

    # load imaginary part of liquid water refractive index and calculate wavelength dependent absorption coefficient
    # __file__ should live at isofit/isofit/inversion/
    
    
    data_dir_path = "../data/"
    path_k = os.path.join(data_dir_path,"k_liquid_water_ice.csv")
    
    #isofit_path = os.path.dirname(os.path.dirname(os.path.dirname(__file__)))
    #path_k = os.path.join(isofit_path, "data", "iop", "k_liquid_water_ice.xlsx")

    # k_wi = pd.read_excel(io=path_k, sheet_name="Sheet1", engine="openpyxl")
    # wl_water, k_water = get_refractive_index(
    #     k_wi=k_wi, a=0, b=982, col_wvl="wvl_6", col_k="T = 20°C"
    # )
    k_wi = pd.read_csv(path_k)
    wl_water, k_water = get_refractive_index(
        k_wi=k_wi, a=0, b=982, col_wvl="wvl_6", col_k="T = 20°C"
    )
    kw = np.interp(x=wl_sel, xp=wl_water, fp=k_water)
    abs_co_w = 4 * np.pi * kw / wl_sel

    rfl_meas_sel = rfl_meas[lw_feature_left : lw_feature_right + 1]

    x_opt = least_squares(
        fun=beer_lambert_model,
        x0=lw_init,
        jac="2-point",
        method="trf",
        bounds=(
            np.array([lw_bounds[ii][0] for ii in range(3)]),
            np.array([lw_bounds[ii][1] for ii in range(3)]),
        ),
        max_nfev=15,
        args=(rfl_meas_sel, wl_sel, abs_co_w),
    )

    solution = x_opt.x

    if return_abs_co:
        return solution, abs_co_w
    else:
        return solution

3.4.1 CWC of a Single Point

Now that we have all of the pieces, we can estimate the CWC of a single pixel using our invert_liquid_water function.

ewt = invert_liquid_water(point.reflectance.values,point.wavelengths.values)
print(f"EWT for ({point.longitude.values:.3f},{point.latitude.values:.3f}): {ewt[0]:.3f} cm")
EWT for (-120.353,34.540): 0.149 cm

3.4.2 CWC of a DataFrame of Points

We can also apply this function to the point data we selected in the previous notebook with the interactive plot.

Read in the csv file we created.

points_df = pd.read_csv("../data/emit_click_data.csv")
points_df
id x y 381.00558 388.4092 395.81583 403.2254 410.638 418.0536 425.47214 ... 2426.4402 2433.8303 2441.2183 2448.6064 2455.9944 2463.3816 2470.7678 2478.153 2485.5386 2492.9238
0 0 -120.569492 34.574552 0.018037 0.017798 0.017567 0.017348 0.017523 0.018065 0.018957 ... 0.034822 0.034534 0.030492 0.031268 0.027023 0.026327 0.024637 0.023704 0.015793 0.013392
1 1 -120.568286 34.591054 0.012029 0.011957 0.011888 0.011826 0.012031 0.012486 0.013151 ... 0.020698 0.022204 0.020115 0.019911 0.016761 0.017437 0.015805 0.016146 0.014353 0.012923
2 2 -120.371752 34.484274 0.020833 0.021323 0.021819 0.022317 0.023065 0.024095 0.025361 ... 0.059581 0.058129 0.053620 0.053622 0.048783 0.044683 0.044064 0.041562 0.034080 0.028254
3 3 -120.181848 34.483303 0.024848 0.028178 0.034840 0.038940 0.043875 0.047232 0.050809 ... 0.131800 0.135400 0.122960 0.118844 0.108201 0.112605 0.106271 0.101222 0.070507 0.059903
4 4 -120.407924 34.591054 0.040804 0.044085 0.048493 0.054733 0.059349 0.064003 0.067749 ... 0.266081 0.263995 0.258724 0.252158 0.228920 0.226200 0.212124 0.181446 0.130708 0.111522
5 5 -120.438851 34.712298 0.015452 0.016034 0.016616 0.017204 0.017848 0.018621 0.019477 ... 0.046427 0.046706 0.043317 0.044288 0.038928 0.037683 0.038662 0.035936 0.034616 0.033542
6 6 -120.512401 34.648230 0.044019 0.039991 0.039511 0.045176 0.045753 0.047011 0.048643 ... 0.081432 0.080438 0.076009 0.076724 0.069166 0.069164 0.066233 0.063523 0.059582 0.056769
7 7 -120.350229 34.626389 0.018255 0.018679 0.019107 0.019550 0.020291 0.021369 0.022694 ... 0.050578 0.047565 0.046092 0.046638 0.039017 0.036911 0.037903 0.033778 0.028663 0.028218
8 8 -120.355052 34.713269 0.017943 0.018290 0.018641 0.019003 0.019586 0.020435 0.021492 ... 0.047776 0.048551 0.045117 0.045340 0.039024 0.038904 0.039077 0.035324 0.030073 0.029010
9 9 -120.243522 34.645318 0.021041 0.021242 0.021448 0.021661 0.022252 0.023214 0.024442 ... 0.049286 0.047242 0.042362 0.045277 0.039692 0.040241 0.037825 0.034632 0.029134 0.029248

10 rows × 288 columns

Now create an array of wavelengths from the column names starting with the first wavelength value (column 3).

# Get wavelength values
wavelength_values = points_df.columns[3::].to_numpy()

Iterate by row through our dataframe, selecting the reflectance values and providing them to the invert_liquid_water function. Afterwards, add an CWC column to our dataframe.

# Creat empty list
ewt_values = []
# Iterate through rows
for _i in points_df.index.to_list():
    # Get reflectance array to pass to function
    rfl_values = points_df.iloc[_i,3::]
    # Use invert liquid water function and append results to list
    ewt_values.append(invert_liquid_water(rfl_values,wavelength_values)[0])    
# Add to our existing dataframe at Column Index 3
points_df.insert(3, "ewt", ewt_values)
points_df
id x y ewt 381.00558 388.4092 395.81583 403.2254 410.638 418.0536 ... 2426.4402 2433.8303 2441.2183 2448.6064 2455.9944 2463.3816 2470.7678 2478.153 2485.5386 2492.9238
0 0 -120.569492 34.574552 2.773249e-01 0.018037 0.017798 0.017567 0.017348 0.017523 0.018065 ... 0.034822 0.034534 0.030492 0.031268 0.027023 0.026327 0.024637 0.023704 0.015793 0.013392
1 1 -120.568286 34.591054 2.495512e-01 0.012029 0.011957 0.011888 0.011826 0.012031 0.012486 ... 0.020698 0.022204 0.020115 0.019911 0.016761 0.017437 0.015805 0.016146 0.014353 0.012923
2 2 -120.371752 34.484274 1.639785e-01 0.020833 0.021323 0.021819 0.022317 0.023065 0.024095 ... 0.059581 0.058129 0.053620 0.053622 0.048783 0.044683 0.044064 0.041562 0.034080 0.028254
3 3 -120.181848 34.483303 5.938623e-02 0.024848 0.028178 0.034840 0.038940 0.043875 0.047232 ... 0.131800 0.135400 0.122960 0.118844 0.108201 0.112605 0.106271 0.101222 0.070507 0.059903
4 4 -120.407924 34.591054 4.001294e-11 0.040804 0.044085 0.048493 0.054733 0.059349 0.064003 ... 0.266081 0.263995 0.258724 0.252158 0.228920 0.226200 0.212124 0.181446 0.130708 0.111522
5 5 -120.438851 34.712298 6.409762e-02 0.015452 0.016034 0.016616 0.017204 0.017848 0.018621 ... 0.046427 0.046706 0.043317 0.044288 0.038928 0.037683 0.038662 0.035936 0.034616 0.033542
6 6 -120.512401 34.648230 2.607293e-10 0.044019 0.039991 0.039511 0.045176 0.045753 0.047011 ... 0.081432 0.080438 0.076009 0.076724 0.069166 0.069164 0.066233 0.063523 0.059582 0.056769
7 7 -120.350229 34.626389 1.505599e-01 0.018255 0.018679 0.019107 0.019550 0.020291 0.021369 ... 0.050578 0.047565 0.046092 0.046638 0.039017 0.036911 0.037903 0.033778 0.028663 0.028218
8 8 -120.355052 34.713269 1.379378e-01 0.017943 0.018290 0.018641 0.019003 0.019586 0.020435 ... 0.047776 0.048551 0.045117 0.045340 0.039024 0.038904 0.039077 0.035324 0.030073 0.029010
9 9 -120.243522 34.645318 1.954748e-01 0.021041 0.021242 0.021448 0.021661 0.022252 0.023214 ... 0.049286 0.047242 0.042362 0.045277 0.039692 0.040241 0.037825 0.034632 0.029134 0.029248

10 rows × 289 columns

For larger sets of point data we would want to do this in parallel. We can also calculate CWC for an area, where we definitly need to utilize parallel processing.

3.5 CWC Calculation of an ROI

In the previous notebook, we subset our region of interest and exported the file. Since the CWC calculation is computationally intensive, it can take a while to process large scenes, so its more efficient to do this spatial subsetting up front. We can use a function included in the ewt_calc.py module to calculate CWC on a cropped image, and create a cloud-optimized GeoTIFF (COG) file containing the results.

Set our input filepaths and output directory.

fp = '../../shared/2023-VITALS-Workshop-AGU/data/cropped/dangermond/rfl/EMIT_L2A_RFL_001_20230401T203751_2309114_002_dangermond.nc'
out_dir = '../data/'
roi_ds = xr.open_dataset('../../shared/2023-VITALS-Workshop-AGU/data/cropped/dangermond/rfl/EMIT_L2A_RFL_001_20230401T203751_2309114_002_dangermond.nc', decode_coords='all')
roi_ds
<xarray.Dataset>
Dimensions:           (wavelengths: 285, latitude: 244, longitude: 262)
Coordinates:
  * wavelengths       (wavelengths) float32 381.0 388.4 ... 2.486e+03 2.493e+03
    fwhm              (wavelengths) float32 ...
    good_wavelengths  (wavelengths) float32 ...
  * latitude          (latitude) float64 34.57 34.57 34.57 ... 34.44 34.44 34.44
  * longitude         (longitude) float64 -120.5 -120.5 -120.5 ... -120.4 -120.4
    elev              (latitude, longitude) float32 ...
    spatial_ref       int64 ...
Data variables:
    reflectance       (latitude, longitude, wavelengths) float32 ...
Attributes: (12/39)
    ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
    summary:                           The Earth Surface Mineral Dust Source ...
    keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
    Conventions:                       CF-1.63
    sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
    instrument:                        EMIT
    ...                                ...
    spatialResolution:                 0.000542232520256367
    spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
    geotransform:                      [-1.20992414e+02  5.42232520e-04 -0.00...
    day_night_flag:                    Day
    title:                             EMIT L2A Estimated Surface Reflectance...
    granule_id:                        EMIT_L2A_RFL_001_20230401T203751_23091...
roi_ds.sel(wavelengths=850,method='nearest').hvplot.image(x='longitude',y='latitude',cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(
    title=f"Dangermond ROI - RFL at 850 nm", xlabel='Longitude',ylabel='Latitude')

Use the calc_ewt function to calculate CWC of the cropped image. This function will also create a COG file containing the CWC results. We can also specify the number of CPUs to use manually with a n_cpu argument, or leave it blank to use all but one of the available CPUs. If we set the return_cwc argument to true, the function will also return the COG.

This will take some time, about 5 minutes, because we’re doing the calculation for roughly 63,000 pixels.

%%time
ds_cwc = calc_ewt(fp, out_dir, return_cwc=True)
2023-12-07 15:11:32,965 WARNING services.py:1996 -- WARNING: The object store is using /tmp instead of /dev/shm because /dev/shm has only 67108864 bytes available. This will harm performance! You may be able to free up space by deleting files in /dev/shm. If you are inside a Docker container, you can increase /dev/shm size by passing '--shm-size=3.05gb' to 'docker run' (or add it to the run_options list in a Ray cluster config). Make sure to set this to more than 30% of available RAM.
2023-12-07 15:11:34,128 INFO worker.py:1664 -- Started a local Ray instance. View the dashboard at 127.0.0.1:8265 
CPU times: user 688 ms, sys: 315 ms, total: 1 s
Wall time: 5min
ds_cwc
<xarray.Dataset>
Dimensions:      (latitude: 244, longitude: 262)
Coordinates:
  * latitude     (latitude) float64 34.57 34.57 34.57 ... 34.44 34.44 34.44
  * longitude    (longitude) float64 -120.5 -120.5 -120.5 ... -120.4 -120.4
    elev         (latitude, longitude) float32 ...
    spatial_ref  int64 ...
Data variables:
    cwc          (latitude, longitude) float64 nan nan nan nan ... nan nan nan
Attributes: (12/13)
    flight_line:            emit20230401t203751_o09114_s000
    time_coverage_start:    2023-04-01T20:37:51+0000
    time_coverage_end:      2023-04-01T20:38:03+0000
    easternmost_longitude:  -119.71545648976226
    northernmost_latitude:  34.7990749308955
    westernmost_longitude:  -120.992414074966
    ...                     ...
    spatialResolution:      0.000542232520256367
    spatial_ref:            GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84...
    geotransform:           [-1.20992414e+02  5.42232520e-04 -0.00000000e+00 ...
    day_night_flag:         Day
    title:                  EMIT Estimated Equivalent Water Thickness (EWT) /...
    granule_id:             EMIT_L2A_RFL_001_20230401T203751_2309114_002

Plot CWC of our ROI.

ds_cwc.hvplot.image(x='longitude',y='latitude',cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(
    title=f"{ds_cwc.cwc.long_name} ({ds_cwc.cwc.units})", xlabel='Longitude',ylabel='Latitude')

Contact Info:

Email: LPDAAC@usgs.gov
Voice: +1-866-573-3222
Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
Website: https://lpdaac.usgs.gov/
Date last modified: 11-28-2023

¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.